Disentangling multiproxy temperature reconstructions from the subtropical North Atlantic

1. Introduction

Many different sea surface temperature (SST) reconstructions based on microfossils and geochemical methods exist for the global ocean, and on time scales from millennia to millions of years. SST reconstructions derived from different proxies that are extracted from the same sediment sample reveal mismatches (e.g., Chapman et al., 1996; Steinke et al., 2008; Jansen et al., 2009; Leduc et al., 2010, 2017; Schneider et al., 2010; Sicre et al., 2013). These differences become evident by comparing, for example, results of the CLIMAP Project Members (1981) and MARGO Project Members (2009) projects that have provided global maps of temperature estimates based on various proxies for the Last Glacial Maximum. While CLIMAP Project Members (1981) reconstructed a glacial warming of the subtropical regions using a combination of different transfer functions, the results of the MARGO Project Members (2009) multiproxy approach may indicate a slight glacial cooling of 2°C for the same region. In addition, multiproxy compilations (Annan and Hargreaves, 2013; Osman et al., 2021) combined with data assimilation partially differ from the MARGO Project Members (2009) dataset.

It is generally argued that discrepancies between different proxy reconstructions are caused by the unique nature of these archives, i.e., the specific ecology of any signal-bearing organism that may be affected by SST changes to different degrees (Jonkers et al., 2010, 2013; Leduc et al., 2010; Schneider et al., 2010; Lohmann et al., 2012; Sicre et al., 2013; Telford et al., 2013; Hessler et al., 2014). Ecological changes may include the seasonal occurrence and varying preferences in dwelling depths of the signal-bearing organisms that are influenced by various oceanographic factors such as vertical and horizontal differences of water mass boundaries, salinity changes, sea ice extent, as well as changes in nutrient availability (Pérez-Folgado et al., 2003; Kisakürek et al., 2008; Storz et al., 2009; Jonkers et al., 2010, 2013; Leduc et al., 2010, 2014; Schneider et al., 2010; Sabbatini et al., 2011; Lohmann et al., 2012; Sicre et al., 2013; Telford et al., 2013; Kretschmer et al., 2016). The effect of variations of these environmental parameters on the marine biota may have a larger effect on the proxy than the temperature change itself (Huybers and Curry, 2006). Thereby, the assumption of a constant relationship between the temperature of ambient seawater and the temperature in given proxy records may be violated (Telford et al., 2013; Laepple and Huybers, 2014).

The discrepancies between SST proxies have been regarded as weaknesses of one or the other method (see summary in MARGO Project Members, 2009). However, careful consideration of ecological differences of the organisms such as planktic foraminifers and coccolithophores may help to understand mismatches between the SST reconstructions. In this study, we reinvestigate a multiproxy temperature dataset, consisting (1) of published, geochemical (Mg/Ca) data from three different planktic foraminifer species (Repschläger et al., 2015, 2017), and (2) U^K′₃₇ temperature reconstructions (Schwab et al., 2012) with (3) published (Weinelt et al., 2015) and new data of faunal assemblages at centennial resolution, and new microfossil based SIMMAX SST reconstructions from the subtropical North Atlantic covering the last 15.5 ka BP. Data are from sediment core GEOFAR KF16/MD08-3180 southeast of the Azores Islands, located at the northern boundary of the subtropical gyre (Figure 1). The coring site is situated between temperate mesotrophic subpolar waters to the north and warmer oligotrophic subtropical waters to the south of the Azores Front (AF).

Phytoplankton production in the North Atlantic increases in spring at low latitudes and successively migrates northward with rising temperatures until late summer (e.g., Henson et al., 2009; Schiebel et al., 2011). Once surface ocean mixing deepens the mixed layer in autumn, biological production increases and a second less pronounced phytoplankton bloom may form (Schiebel et al., 2001). The strength, timing, and location of phytoplankton production in Azores waters is highly sensitive to changes in the strength of the North Atlantic oscillation (NAO) (Henson et al., 2009). Spring phytoplankton production starts earlier during negative than during positive NAO phases. This strong coupling of the phytoplankton bloom to the NAO indicates the strong coupling between the AF region and the North Atlantic climate on short timescales. On longer timescales, the Acores Front Current System is assumed to affect the stability of Atlantic Meridional Overturning Circulation (AMOC) as it acts as a potential reservoir for heat during glacial and deglacial weak phases of the AMOC (e.g., Repschläger et al., 2015, 2017).

On glacial to interglacial time scales, the position of the AF is thought to have changed latitude (Schiebel et al., 2002a). Other records from the Azores coring site indicate surface water cooling and southward movement of the AF during cold Heinrich event 1 (H1) and Younger Dryas, YD (Schwab et al., 2012; Repschläger et al., 2015). A relatively low amplitude in the U^K′₃₇ temperature record compared to the Iberian Margin is observed for the deglacial period and explained with a shift of the main alkenone production from spring to summer (Schwab et al., 2012) during cold periods. Surface and subsurface temperature reconstructions from different foraminifer species have been explained by heat storage in subsurface waters (Repschläger et al., 2015, 2017). These results indicate that the ecology of the different proxy bearing organisms at the coring site affects the observed temperature signal.

In the following,

(1) We test the behavior of different temperature proxies over the last deglaciation under consideration of archive-specific differences in ecology. To better explain the discrepancies between individual SST reconstructions, and to gain better insight into the deglacial and Holocene changes in surface hydrology of the Azores Current (AC) with respect to (a) seasonal temperature changes, (b) surface ocean stratification, and (c) water mass changes.

(2) We assess the plausibility of the SST reconstructions comparing our SST with published SST reconstructions (Obbink et al., 2010; Thornalley et al., 2010, 2011; Repschläger et al., 2015) and fresh water-forced simulation (i.e., salinity) experiments.

(3) To assess local salinity changes, classically, ice volume and temperature effects are subtracted from stable oxygen isotope measurements performed on surface-dwelling foraminifers (e.g., Shackleton, 1974; Duplessy et al., 1991; Waelbroeck et al., 2002). For this purpose, temperature reconstructions from Mg/Ca, U^K′₃₇, and SIMMAX are used and combined with the δ¹⁸O signal of surface and subsurface dwelling planktic foraminifer species of varying phenology (Schiebel and Hemleben, 2000). Under the assumption that the different temperature reconstructions are potentially reflecting different seasons and water depths, the combination of temperature and δ¹⁸O values from different proxy bearing organisms might have a strong impact on the calculations and may overprint the actual calculated δ¹⁸O_WATER (δ¹⁸O_w) values. This effect will be tested in the last paragraph of this paper.

2. Regional setting

The surface hydrography at the GEOFAR KF16/MD08-3180 coring site in the subtropical North Atlantic is governed by the Azores Front-Current System (Figure 1). The modern Acores Current (AC) branches off from the North Atlantic Current (NAC) between 30 and 40°N and meanders east, inducing mesoscale eddies (e.g., Käse and Siedler, 1982). Within this current system, the AF separates oligotrophic Subtropical Mode Water (SMW), and mesotrophic Eastern North Atlantic Central Water (ENACW) (Käse and Siedler, 1982; McCartney and Talley, 1982; Mason and Oliveira, 2006; Liu and Tanhua, 2019a,b). In newer publications, these water masses are also named tropical Eastern North Atlantic Central Water (ENACW_T) and subpolar North Atlantic Central Water (ENACW_P) (e.g., Mason and Oliveira, 2006; Liu and Tanhua, 2019a,b). Both SMW and ENACW are surface water masses (e.g., Liu and Tanhua, 2019b). The AF is defined by the position of the southward dipping 15°C-isotherm (Gould, 1985). Being a subsurface front overlain by the surface waters with relatively homogeneous SST (Alves et al., 2002), changes in SST across the AF cannot be used to trace the position of the AF. The divide between different surface dwelling planktic foraminifer assemblages is located between the SMW and ENACW to the north of the AF (Ottens, 1991; Schiebel et al., 2002a; Rogerson et al., 2004) at about 34–38°N.

Under modern conditions, main phytoplankton production in the vicinity of the Azores occurs at increasing light-levels in early spring (March) and moves north to the transitional to subpolar NA (40 and 60°N) following photic conditions until June (Obata et al., 1996; Longhurst, 2007; Henson et al., 2009). Coccolithophore production in the mid-latitude North Atlantic may persist beyond maximum production of other phytoplankton groups, and enhanced coccolithophore diversity (i.e., number of species) may characterize the late spring to early summer succession in surface waters (e.g., Haidar and Thierstein, 2001; Schiebel et al., 2011; Hopkins et al., 2015; Poulton et al., 2017).

3. Materials and methods

A high-resolution deglacial record of the combined cores GEOFAR KF16 (37°59.59′N 31°7.698′W, 3050 m water depth) and MD08-3180 (37°59.99′N 31°8.07′W, 3064 m water depth) from a coring site situated to the southwest of the Azores at the northern boundary of the AF has been analyzed (Richter, 1998; Kissel et al., 2008; Schwab et al., 2012; Repschläger et al., 2015). A previously published age model is based on 15 Accelerator Mass Spectrometer (AMS) radiocarbon ages of monospecific planktic foraminifer samples, and refined by tuning the alkenone SST record to NGRIP ice core data (Schwab et al., 2012). The original age model has been modified by using the ¹⁴C AMS plateau tuning technique (Sarnthein et al., 2015). Planktic foraminifer tests are well preserved as shown by the occurrence of pteropods during the YD. Planktic foraminifer assemblage data were previously published by Repschläger et al. (2015) and Weinelt et al. (2015), and supplemented with new Holocene data (97 samples) in centennial resolution (see Supplementary Table 1 for details). All planktic foraminifer counts were performed on a split of at least 300 individuals from the size fraction >150 μm following the taxonomy of Pflaumann et al. (1996), which is based on Bé and Tolderlund (1971).

Temperature reconstructions are compared for three different methods (1) planktic foraminifer transfer functions (SIMMAX) (Pflaumann et al., 1996), (2) alkenone SST reconstructions (U^K′₃₇ SST) (Schwab et al., 2013), and (3) Mg/Ca temperature reconstructions (Repschläger et al., 2015, 2017). For the Holocene, Mg/Ca and SIMMAX data were obtained from the same samples taken from core GEOFAR K16 and U^K′₃₇ records stem from core MD08-3180. For Bølling-Allerød (BA) and YD, all analyses were performed on samples from core MD08-3180. Both Mg/Ca and stable isotope analyses was performed on specimens of the size range >350 μm.

SIMMAX (Pflaumann et al., 1996) were calculated using published planktic foraminifer assemblage data (Repschläger et al., 2015; Weinelt et al., 2015) and supplemented with unpublished Holocene planktic foraminifer assemblage data at decadal to centennial resolution. The used reference database SIMMAX28 consists of 947 core top samples that are linked to SST from the Levitus94 database (Levitus and Boyer, 1994). To calculate distant and non-distant weighted summer and winter SSTs the best ten analogs have been used. The minimum correlation coefficient of the default SIMMAX setting was reduced from 0.9 (90%) to 0.0 (0%) as the use of default settings resulted in non-analog situations for more than 30% of the samples. For the adjusted setting, the correlation coefficient is >0.92 for all Holocene samples except samples around 8.2 ka BP, 9.3–9.6 and 10.4 ka BP with lower similarities. Samples older than 11.4 ka BP show correlations between 0.7 and 0.9. Within this study, not distant weighted results were used, as core top values better reproduce the modern situation. The difference to the distant weighted method does not exceed 0.5°C. The standard deviation of the SIMMAX SST reconstruction is ± 1°C.

The alkenone unsaturation index U^K′₃₇ (Prahl and Wakeham, 1987) is calibrated against temperature using the calibration of Muñller et al. (1998) and was published by Schwab et al. (2013). The Mg/Ca data from the size fraction >315 μm were previously published by Repschläger et al. (2015, 2017), the standard deviation of SST calculation for this core is ± 0.2°C (Schwab et al., 2013) obtained from replicate measurements. The systematic error of the calibration used in this study is ± 1°C. Mg/Ca based temperature estimates from surface (Globigerinoides ruber, Globigerina bulloides) and subsurface dwelling (Globorotalia truncatulinoides) planktic foraminifers are calibrated with the equations given by Anand et al. (2003) (general equation for mixed planktic foraminifera species, eq. 5) and Cléroux et al. (2008) (eq. 1). These calibrations have been chosen from a variety of existing calibrations as their results give the best match for Late Holocene (LH) temperatures with the modern conditions at the AF (see Supplementary Table 1). Duplicate measurements of downcore samples result in a combined measurement and calibration with a standard deviation of ± 1°C for G. ruber and G. bulloides, and ± 2°C for G. truncatulinoides. Mg/Ca measurements on G. ruber (Mg/Ca SST_Gr) have only been possible for the warm Allerød and Holocene. Low numbers of tests of G. ruber impede analyses of the cold intervals YD and H1. Mg/Ca was measured on G. bulloides over the late BA to late YD time interval. During the BA and Intra Allerød Cold Period (IACP) time intervals, several results had to be discarded due to mineral coatings indicated by anomalously high Fe content. Full analytical methods and errors of temperature reconstruction are presented in the original publications (Schwab et al., 2013; Repschläger et al., 2015, 2017; Weinelt et al., 2015). For transparency, all data compiled for this study can be found with original citations in the Supplementary Tables 1–6.

Abundances of planktic foraminifer species are used as indicators for additional evidence on hydrological changes at the Azores coring site: G. ruber [all white morphotypes combined for their common ecological significance (Yamasaki et al., 2022)] are tropical to subtropical surface- dwellers that typically thrive in waters warmer than 18°C (Bradshaw, 1957). Samples from core tops (Pflaumann et al., 1996), sediment traps, and plankton tows (Schiebel et al., 2002b) indicate a dominance of its abundance during late spring at 32°N, and a change to late summer conditions further to the North, and are associated with oligotrophic late summer conditions (Storz et al., 2009). Here, we assume G. ruber a is an indicator for the STG waters, with its variability in relative abundance reflecting changes in the AF position. Additionally, we propose that the Mg/Ca record of G. ruber reflects summer temperatures.

Deep-dwelling G. truncatulinoides are most abundant south of the AC system within the SMW realm. G. truncatulinoides is assumed to have an annual life cycle and to add a secondary calcite at thermocline depth before ascending for reproduction in the Azores region (Hemleben et al., 1989; Schiebel et al., 2002b). Surface and subsurface dwelling Neogloboquadrina incompta and Globorotalia scitula (Hemleben et al., 1989; Ottens, 1991; Schiebel et al., 2002a,b; Rebotim et al., 2017, 2019) typically occur north of the AF within the ENACW. Enhanced relative abundances of G. scitula are indicative of the northern water bodies (Schiebel et al., 2002a).

Surface to subsurface dwelling Turborotalita quinqueloba and Neogloboquadrina pachyderma (Rebotim et al., 2017, 2019) are most abundant in the Nordic Seas and at the Polar Front being valid indicators for subpolar to polar conditions (Bé and Tolderlund, 1971; Ottens and Nederbragt, 1992; Pflaumann et al., 1996; Simstich et al., 2003). T. quinqueloba is also associated to hydrologic fronts and occurs in high numbers at the warm side of the Arctic front (Johannessen et al., 1994; Husum and Hald, 2012). The cold-water species N. pachyderma (previously N. pachyderma left coiling) indicates the influence of subpolar waters at the coring site.

Globigerina bulloides and Globigerinita glutinata are most abundant in subpolar to transitional waters showing a bimodal abundance peaking in late spring and autumn (Schiebel and Hemleben, 2000; Schiebel et al., 2001; Jonkers and Kučera, 2015).

The planktic foraminifer δ¹⁸O_c record (Repschläger et al., 2015, 2017; Balmer and Sarnthein, 2018) is corrected for the temperature effect, using SIMMAX, Mg/Ca and U^K′₃₇ temperature reconstructions and the Shackleton (1974) equation. This results in estimates for changes in past δ¹⁸O_w composition. This signal has been further corrected for the ice-volume effect, using the relative sea level composite curve of Austermann et al. (2013) to calculate δ¹⁸O_w^–ice.

4. Results

The combined spliced sediment record spans the last 15.5 ka BP with a multidecadal to centennial resolution. It includes the last deglaciation with warm periods such as the BA and the Holocene, and cold periods such as H1 and YD events. The short-term cold Greenland Interstadial 1b at 13.1 ka BP, also known as Intra Allerød Cold Period (IACP) (Rasmussen et al., 2006, and references therein), and the 8.2-kyr-BP-Event interrupt the warm BA and Holocene, respectively. According to the SIMMAX and U^K′₃₇ derived SST estimates, the general temperature pattern follows the known deglacial succession with warm BA and Holocene temperatures, and cold H1 and YD temperatures (Bard, 2001a) in the North Atlantic. Yet, distinct differences are observed between records (Figure 2). In general, SIMMAX summer and winter temperatures show the same pattern of change with 5°C colder temperatures during winter in average. Smaller winter than summer cooling amplitudes during cold H1 and late YD/early Holocene are observed.

4.1. Heinrich event 1 (H1)

Coldest surface temperatures are observed during H1 in all surface temperature records. SIMMAX summer SST (SIMMAX SST_su) (Figure 2A) drop to 9°C, and thus indicate a summer cooling of 13°C with respect to modern SSTs. U^K′₃₇ SSTs indicate temperatures around 12.5°C (Figure 2B), with a cooling amplitude of 7.5°C that is less pronounced than in the SIMMAX record. Although calibrated to reflect annual mean SSTs, U^K′₃₇ temperatures during H1 are warmer than planktic foraminifer-based summer SSTs.

Globorotalia truncatulinoides based subsurface Mg/Ca temperatures (Mg/Ca T_Gtrunc) indicate temperatures of 14°C for H1, Figures 2C, 3E, H. Numbers of subtropical planktic foraminifers are low during H1, with relative abundances of G. ruber and G. truncatulinoides at <2% (Figures 3A, B), while the relative abundances of polar and subpolar species such as N. pachyderma (5%), T. quinqueloba (15–20%), and N. incompta (15–30%) (Figures 3C–G) are higher and increase toward the end of H1. Also, relative abundances of G. scitula are elevated (10–20%) in H1 in comparison to the preceding and following time-intervals (Figure 3C).

4.2. Bølling-Allerød (BA)

Bølling-Allerød interval is 5 to 9°C warmer than the H1 event with average U^K′₃₇ and SIMMAX SST_su of 19 to 20°C and match modern annual mean SST of 19°C (Figures 2A, B). SIMMAX SST_wi remain 5°C colder. Mg/Ca SST_Gr are only available for the late Allerød, indicating SST of 21.5°C close to LH SST. G. bulloides Mg/Ca SST (Mg/Ca SST_Gb) are similar to U^K′₃₇ derived SSTs, reaching modern winter SSTs of 17°C from Bølling to early Allerød (14.4–14.2 ka) (Figure 2C). Here, the IACP interval in which Mg/Ca SST_Gb are overprinted by mineral coatings are excluded. Subsurface Mg/Ca temperatures of G. truncatulinoides (Mg/Ca T_Gtrunc) decrease from 16°C at the beginning of BA to 14°C at the onset of IACP (Figure 2C). Relative abundances of G. ruber were at intermediate values (10%) during BA (Figure 3A). The relative abundances of N. incompta (Figure 3F) show a pronounced peak (30–35%) at the onset of BA (12.6 to 12.4 ka BP). Relative abundances of polar/subpolar and temperate species N. pachyderma (0–2%) (Figure 3H), T. quinqueloba (5%) (Figure 3G), and G. glutinata (5–7%) (Figure 3F) are low. Relative abundances of G. scitula decrease from about 8 to 2% in early BA, at 14.5 ka BP, and remain low until IACP (Figure 3C).

The IACP, centered at 13.1 ka BP, interrupts warm BA with a pronounced cooling in the SIMMAX SST_su to 10°C, i.e., temperatures close to those estimated for the H1 event. The strong cooling in SIMMAX SST estimates coincides with a cooling signal of 1.5°C in the U^K′₃₇ record from 18.5 to 17°C, which is represented by a single data point only. Relative abundances of T. quinqueloba and G. scitula sharply increase to 30 and 14%, respectively (Figures 3C, G). After IACP, SSTs and abundance records return to average BA values, with the exception of relative abundances of G. ruber a that remain low at <2% until the end of YD.

4.3. Younger Dryas (YD)

At the beginning of YD, 4°C cooling is observed in the U^K′₃₇ SST and Mg/Ca SST_Gb records (Figures 2B, C). After the initial cooling, both records range between 13 and 15°C. and Mg/Ca T_Gtrunc remain relatively warm with average values of about 14°C. Mg/Ca SST_Gb, U^K′₃₇ SST, and Mg/Ca T_Gtrunc provide similar temperatures that match modern winter temperatures and imply a collapse of the thermal surface water stratification. The exception in this pattern is the SST_su record that remains relatively warm (16°C) until 12.2 ka and cools down by 2 to 14°C thereafter (Figure 2A). Over the entire YD, SST_wi are 4°C colder than SST_su. During late YD, SIMMAX SST_su cools down by 4 to 14°C but still remains 2°C warmer than the U^K′₃₇ SST. The relative abundances of G. ruber a, G. truncatulinoides, N. pachyderma, and G. scitula are very low (<2%) during YD. The relative abundances of N. incompta show two peaks, one at the onset of the YD (37%) and one (40%) in late YD (Figure 3). Relative abundances of T. quinqueloba (Figure 3G) are slightly decreasing from 15 to 10%, and abundances of G. bulloides fluctuate by 15 around 30%. Abundances of G. glutinata vary between about 25 and 10% (Figure 3E).

4.4. Early Holocene

The general warming trend into the Early Holocene from average temperatures of 14°C during YD to 23°C at 9.5 ka BP in the SIMMAX SST is characterized by pronounced and rapid decadal fluctuations with magnitudes of about 4°C (Figure 2A). This short-term climate variability is not visible in the U^K′₃₇ SST record (Figure 2B). U^K′₃₇ SST rise in the beginning of the Early Holocene and reach modern annual mean SST of 19°C at 11.5 ka BP, and stay constant until the 8.2-ka-BP-Event. Mg/Ca SST_Gr indicates relatively stable SST of 23°C from 10.5 ka BP on. Short-term variability cannot be tested in the Mg/Ca SST_Gr and Mg/Ca SST_Gb records due to missing data for these events (Figures 2A, C). The Mg/Ca T_Gtrunc cool down to 12°C, reaching H1 values at the beginning of Early Holocene and show a long-term gradual increase into the 8.2-ka-BP-Event to 17°C. Relative abundances of G. ruber are at intermediate values (10%). Relative abundances of G. truncatulinoides fluctuate between 1 and 4% and are slightly increased compared to the YD to H1 interval before (Figure 3). Relative abundances of N. pachyderma (<2%) and G. scitula (1–6%) are low, N. incompta relative abundances steadily decrease from 40% at 11.6 ka and reach Holocene mean values of about 10% at 10 ka BP (Figures 3C, F, H). Relative abundances of T. quinqueloba (Figure 3G) slightly decrease from 10 to 5% and G. bulloides relative abundances vary by 10% around a mean value of 30% (Figure 3D). G. glutinata (Figure 3E) show similar variability of ± 9% around a mean of 15% with an inverse distribution to G. bulloides.

4.5. Mid- to Late Holocene

From the Mid to Late Holocene, SIMMAX SST_su/wi are close to modern summer (23°C) and winter surface temperatures (16°C). Over the past 9 ka, temperatures vary by 1°C being slightly warmer from 6.5 and 4 ka BP SST_su/wi and cooler from 9 to 6.5 ka BP and 4 to 0.5 ka BP. The temperature difference between summer and winter amounts to about 7°C during the entire Holocene (Figure 2A). The U^K′₃₇ SST record resembles the annual mean SSTs over the Holocene (19°C) with an SST maximum between 11 and 8 ka BP, following maximum northern hemisphere summer insolation (Figures 2B, 4). Mg/Ca SST_Gr indicate SSTs of 21–22°C, slightly below modern summer SSTs (Figure 2C). Mg/Ca T_Gtrunc fluctuates around an average of 14°C coinciding with modern subsurface temperatures at 200 m water depth.

G. ruber reach maximum relative abundances of about 25% between 8 and 4 ka BP (Figure 3A). After 4 ka BP, the abundances stabilize on LH values between 14 and 18%. The relative abundances of G. truncatulinoides, G. scitula, and N. incompta increase from 7 ka (2%) until the end of the Holocene (10%) (Figures 3B, C, F). Abundances of G. bulloides strongly fluctuate around an average of 32% (Figure 3D), while relative abundances of the polar and subpolar species N. pachyderma (<2%) and T. quinqueloba (<5%) are low during the entire Mid- to Late Holocene (Figures 3G, H). Abundances of G. glutinata fluctuate around 7% and seem to follow a systematics different from G. bulloides (Figures 3D, E).

5. Discussion

5.1. Reliability of records

Mg/Ca temperature reconstructions from planktic foraminifer shells are affected by seasonal changes in productivity, species specific water depth habitats, and the analyzed test size fractions (Skinner and Elderfield, 2005; Chapman, 2010; Schneider et al., 2010; Lombard et al., 2011; Friedrich et al., 2012; Lohmann et al., 2012; Rebotim et al., 2017), as well as changes in the chemistry of ambient seawater including salinity (Kisakürek et al., 2008; Mathien-Blard and Bassinot, 2009; Sabbatini et al., 2011), and carbonate ion concentration (Arbuszewski et al., 2010; Evans et al., 2016; Gray et al., 2018). Furthermore, trace element coatings of planktic foraminifer shells can strongly affect the Mg/Ca signature (Barker et al., 2003; Bian and Martin, 2010). Reduced reliability of the Mg/Ca temperature calibration at the cold end of the calibration dataset might lead to an overestimation of the SSTs lower than about 5°C (Meland et al., 2005; Kozdon et al., 2009). These factors will be addressed in the following.

A variety of calibration functions to transfer Mg/Ca from foraminifer shell calcite into temperature are available. These functions are based on cultured individuals (Nürnberg et al., 1996; Kisakürek et al., 2008), sediment trap data (e.g., Anand et al., 2003; McConnell and Thunell, 2005), and core top data (e.g., Elderfield and Ganssen, 2000; Lea et al., 2000; Dekens et al., 2002; Cléroux et al., 2007, 2008; Regenberg et al., 2007) including multispecies (Elderfield and Ganssen, 2000; Anand et al., 2003) and monospecific approaches (Elderfield and Ganssen, 2000; Lea et al., 2000; Dekens et al., 2002; Cléroux et al., 2007, 2008). The use of different calibration functions results in offsets between the calculated Mg/Ca temperatures (Supplementary Figure 1). Thereby, differences in calibration functions affect the slope and intercept of the calibrations, leading to differences in the absolute temperatures of 1–4°C for G. ruber, 0.1 –1.8°C for G. bulloides, and 2.8–7.9°C for G. truncatulinoides (Supplementary Figure 1). Within this range of temperature estimates, some result in unreasonably high or low temperature reconstructions for the core top samples at our coring site, and can thus be excluded. The differences of absolute temperature estimates of residual calibrations are within the range of the calibration uncertainty. From all calibrations, downcore amplitudes of change in Mg/Ca temperatures from G. ruber and G. bulloides differ by less than 1°C, and by 0.5–1.4°C in G. truncatulinoides. These differences in amplitude are well within the standard deviation of the calibrations.

The differences in absolute temperatures have been explained by various effects. Differences in ontogenetic stage are shown to have a major effect on the Mg/Ca and δ¹⁸O values in planktic foraminifers (Friedrich et al., 2012). This is especially important for shell sizes below 250 μm. The shell sizes used in this study are >315 μm and calibration functions obtained from similar shell sizes were used. Thus, shell size dependant offsets can be regarded as minor and will not be addressed in the following discussion.

Differences between core top, culture, and sediment trap calibrations are explained by early diagenesis in the sediments (Regenberg et al., 2007; Hertzberg and Schmidt, 2013). By addition of gametogenic calcite on the shells that only becomes evident in foraminifers that underwent the full life cycle. These differences may not occur in plankton tow samples (Jentzen et al., 2018). In addition, many core top calibrations infer the dwelling depth of foraminifers by comparing foraminifera δ¹⁸O values with calculated δ¹⁸O calcite (δ¹⁸O_c) data (e.g., Regenberg et al., 2009). These values are calculated in two steps. Due to a lack of δ¹⁸O measurements in the water column (δ¹⁸O_w), δ¹⁸O_w is calculated from salinity measurements using local salinity- δ¹⁸O_w conversions. In a second step, δ¹⁸O_c values are calculated using species specific equations that convert δ¹⁸O_w under consideration of water temperature into theoretical δ¹⁸O_c for all water depths. The comparison of δ¹⁸O_c with the actual δ¹⁸O measurements in the foraminifers is used to infer the dwelling depth of the foraminifer. This practise has several weaknesses: The conversion from salinity into theoretical δ¹⁸O_w data neglects seasonal variability, and water mass differences in the δ¹⁸O_w to salinity relationships (McConnell et al., 2009). For conversion of δ¹⁸O_w into δ¹⁸O_c, a wide range of different formula are available, resulting, for example, in calculated calcification depths of the upper 80 m for G. ruber and 200 m for G. truncatulinoides (Regenberg et al., 2009). Furthermore, the use of annual data for the calibration neglects differences in the seasonal abundance of foraminifers, and associated changes in ambient water temperature and salinity.

Differences in calculated calcification depths transfer, for example, to 4°C temperature differences in the subtropical North Atlantic that are reflected in the absolute SST difference calculated by the various Mg/Ca calibrations. The effect of (neglected) seasonality in calibration functions may account for temperature differences of up to 4°C for G. ruber and 2°C for G. truncatulinoides (Regenberg et al., 2009). Furthermore, most calibrations are focused on ocean basins rather than using global datasets to account for differences in ocean chemistry (e.g., Lea et al., 2000).

Given the multiple effects that lead to the observed differences between sediment traps, living foraminifers, and core top samples, as well as differences between ocean basins, in the following we only consider calibrations from the Atlantic Ocean basin. In agreement with other studies, we chose the calibration function of Anand et al. (2003) that represents the closest match between the modern observations and core top data considering that G. ruber reflects late summer conditions in the upper 30–50 m (Schiebel et al., 2002b; Storz et al., 2009; Rebotim et al., 2017). The evaluation of Mg/Ca calibrations on G. bulloides are impeded by the lack of core top data. In this case, we chose the calibration from the same North Atlantic core top dataset as for G. ruber (Anand et al., 2003). Based on the plankton tow and sediment trap data we assume that G. truncatulinoides match the calibration of Cléroux et al. (2007).

In addition to early diagenesis, an underestimated influence of salinity and carbonate chemistry (pH and CO₃^2– concentration) on the Mg/Ca incorporation (Gray et al., 2018; Gray and Evans, 2019) might have an effect on the SST reconstructions. These effects have to be addressed for the downcore samples, though core top and Late Holocene SST reconstructions match well with modern conditions. However, SEM images of the coccolithophores (Schwab, 2013) show good carbonate preservation at our coring site MD08-3180, and argue against effects of early diagenesis (Schwab, 2013). Following Hertzberg and Schmidt (2013), we assume that carbonate dissolution had no major effect on the preservation of foraminifer shells at our coring site as well.

At the coring site southwest of the Azores, glacial-to-interglacial variability in salinity and pH may account for decreased surface water temperatures of 1–1.5°C for G. ruber and G. bulloides, in comparison to model data shown by Gray and Evans (2019). This discrepancy is slightly exceeding the already considered combined error of measurement and calibration. Corrections for salinity and pH are hampered by the lack of independant salinity estimates, and pH reconstructions bear the risk of additional uncertainty arising from uncertainties in the pH reconstructions of ± 0.1 pH units (Raitzsch et al., 2018). The latter is in the same range as the uncertainty when using uncorrected data. Therefore, we did not apply any salinity and pH correction to our data. However, underestimation of the effects of pH and salinity could result in an underestimation Mg/Ca SST of 1°C, and would lead to an overestimation of the glacial-to-interglacial SST amplitude in all Mg/Ca SST estimates, and underestimate the difference to the alkenone SST record. To conclude, we assume that interspecies Mg/Ca SST temperature amplitudes are robust, whereas the error in absolute SST estimates might be higher than the ± 1°C standard deviation of Mg/Ca measurements in G. ruber and G. bulloides, and differences to the alkenone SST record might be underestimated.

To transfer alkenone records into SST, two different indices can be used: The U^K₃₇ index including the 37:2, 37:3, 37:4 alkenones, while the U^K′₃₇ index includes 37:2 and 37:3 only. Ocean-basin specific calibrations are available. The effect of 37:4 on our SST calculations has already been tested by Schwab et al. (2012). The U^K₃₇ Index based SST calculations, as well as the North Atlantic calibration set (Muñller et al., 1998) in H1 provide 1°C colder SSTs than U^K′₃₇ based SST reconstructions. These calibration differences are within the standard deviation of ± 1°C and can be regarded as a minor effect.

Although calibrated to defined water depths (here 10 m) in the surface mixed layer, SST reconstructions with the modern analog technique using SIMMAX (Pflaumann et al., 1996) might be affected by a strong imprint of subsurface (below thermocline) water mass conditions (Adloff et al., 2011; Telford et al., 2013). This is caused by the addition of species from different water depths in fossil planktic foraminifer assemblages, thus integrating the thermal structure of the whole upper water column. As long as the thermal structure of the upper ocean is stable, planktic foraminifer assemblages reflect SST. A change in the structure of the water column may have a stronger effect on the subsurface dwelling species than surface dwellers, and lead to a bias in the resulting temperature signals toward a subsurface water signal. This bias is especially important in regions with pronounced oceanographic changes between glacial and interglacial times that lack analog situations in the calibration datasets (Storz et al., 2009; Adloff et al., 2011; Telford et al., 2013). One of these regions is the northern boundary of the North Atlantic Subtropical Gyre (STG) where strong fluctuations of the AF are observed (Repschläger et al., 2015) that cause non-analog situations.

Biases in the transfer functions by non-analog conditions to the calibration dataset are highly possible given the non-analog character of our samples (see section “3. Materials and methods”). Thus, in the following, we rather address the feasibility of calculated SST estimates than discuss absolute SSTs. Estimates of winter temperatures and seasonal differences are most probably not independent due to a strong dependence on the summer dataset of calibration (Kucera et al., 2005). Thus, we assume that winter temperatures are less reliable. The differences in seasonal temperature amplitudes obtained from SIMMAX as well as SIMMAX winter SST estimates are not discussed in the following, and are presented for completeness of the classical SIMMAX approach only (Figures 2, 6).

5.2. Differences between Mg/Ca SST, SIMMAX SST, and U^K′₃₇ temperatures during the warm Late Holocene and BA time-intervals

Reconstructed average Mg/Ca SST_Gr of 20.8°C during BA, and 21.8°C during Late Holocene differ in 1°C, and can be assumed similar, when considering an error of ± 1°C in SST estimates (Supplementary Figure 3).

Average U^K′₃₇ SST of 18.5°C in BA are similar to Mg/Ca SST_Gb, yet 2.3°C cooler than Mg/Ca SST_Gr. Although both Mg/Ca SST_Gr and U^K′₃₇ SST warm to 21.8/21.9 and to 19.2/18.9°C, respectively, during early/late Holocene, the difference between average Mg/Ca SST_Gr and U^K′₃₇ SST reconstructions increases to 2.7°C in early Holocene and to about 2.9°C during the Late Holocene. Average Mg/Ca SST_Gr are 3.3°C warmer than SIMMAX SST_su during BA, and are similar during Late Holocene. To explain these differences, we need to understand the differences in ecology of the proxy bearing organisms.

Observational data show maximum production of G. ruber in the modern subtropical to temperate North Atlantic following overall enhanced primary productivity in late spring, and maximum relative abundances in late summer (Hemleben et al., 1989; Wolfteich, 1994; Schiebel et al., 2002b; Storz et al., 2009; Salmon et al., 2015; Schiebel and Hemleben, 2017). This finding is corroborated by results from ecological planktic foraminifer models for the modern situation (Kretschmer et al., 2018), and indicate similar phenology of G. ruber during glacial and interglacial times (Fraile et al., 2009; Lombard et al., 2011). Bearing photo-symbiosis, the dwelling depth of G. ruber is limited to the upper photic zone (e.g., Hemleben et al., 1989; Rebotim et al., 2017) with an average dwelling depth of 50 m in the Azores region (Schiebel et al., 2002a,b). Mg/Ca SST_Gr reconstructions of 21.5°C for late Holocene match well with late July SST at the coring site. Relative stable Late Holocene and BA Mg/Ca SST_Gr of 21.5 to 21.0°C indicate only small changes in water temperatures recorded by G. ruber for both time-intervals. This can either be associated to similar SST conditions, or to a change in the abundance maximum of G. ruber toward late August in BA. To test both possibilities, we compare the Mg/Ca SST_Gr record to the other SST calculations (Figures 4, 6).

The 2.8°C warmer Mg/Ca SST_Gr than U^K′₃₇ SSTs during the BA and Early Holocene are readily explained by the differential ecologies of the foraminifer species G. ruber and alkenone producing coccolithophores, respectively (Figure 4 and Supplementary Figure 3). Main coccolithophore production in the subtropical to subpolar North Atlantic occurs slightly after the observed annual chlorophyll maximum at the end of March until early June (Henson et al., 2009; Narciso et al., 2016). Coccolithophores need light to perform photosynthesis and thus are bound to the euphotic zone. U^K′₃₇ SST reconstructions match well with modern SST during end of May at 0–50 m water depth, Figure 1B, i.e., the average depth habitat of the regionally most abundant coccolithophore species Emiliania huxleyi (Haidar and Thierstein, 2001; Schiebel et al., 2011). Thus, we assume that U^K′₃₇ SST represent late May surface temperatures in the Azores region. While U^K′₃₇ is calibrated to match mean annual SST, it displays the spring SST in the Azores region. The difference between late May and July/August SST at the Azores coring site explains the general difference in Mg/Ca SST_Gr and U^K′₃₇ SST of 2.8°C. The decrease in the SST difference during BA to 2.3°C most probably is related to a slight shift of the phenology of the coccolithophore bloom toward later in the year or to slightly warmer spring temperatures during BA. A phenology change is proposed to occur during positive phases of NAO (Henson et al., 2009) and may argue for a NAO⁺ like pattern during BA that needs to be tested in future projects.

Mg/Ca SST_Gb match well with U^K′₃₇ SSTs during BA and YD. Observations and modeling studies (Schiebel et al., 1997, 2002a; Rebotim et al., 2017; Kretschmer et al., 2018) suggest that G. bulloides mainly dwells in the upper 60 to 100 m in the temperate North Atlantic. As G. bulloides does not form major secondary calcite, it can be assumed that calcification depth does not differ from this water depth interval (Schiebel et al., 1997). Average Mg/Ca SST_Gb of 15.5°C in BA match well with modern SST observed in the upper 60 to 100 m water depth at the end of April. The depth habitat of G. bulloides is mainly controlled by food availability that coincides with the distribution of phytoplankton (Schiebel et al., 2001; Wilke et al., 2009). G. bulloides is feeding on diatoms from the early phytoplankton spring bloom that occurs slightly before major mass production of coccolithophores at the end of spring (Schiebel et al., 1995; Poulton et al., 2017). Thus, we assume that Mg/Ca SST_Gb and U^K′₃₇ SSTs reflect the same season, i.e., spring, with a short time lag between the maximum abundances of both species G. bulloides and E. huxleyi.

Comparison of Mg/Ca SST_Gb and U^K′₃₇ SSTs with the Mg/Ca T_Gtrunc shows a unique temperature development between surface and subsurface waters. Mg/Ca T_Gtrunc during BA is warmer than the Mg/Ca SST_Gb and reaches similar values as under modern conditions at 100 m water depth (Repschläger et al., 2015). In combination with stable isotope data this phenomenon has been ascribed to subsurface warming, which may also have had an influence on temperatures of the overlying water column and may explain the reduced difference between Mg/Ca SST_Gr and U^K′₃₇ SSTs during BA in comparison to the Holocene.

Despite small temperature differences between Late Holocene and BA reconstructed from Mg/Ca (1°C) and U^K′₃₇ (0.4°C), the SIMMAX SST record indicates a 2.5°C difference between both time intervals (Supplementary Figure 3). During Mid and Late Holocene, SIMMAX SST_su (Supplementary Figure 3) match well with Mg/Ca SST_Gr. In contrast SIMMAX SST_su are 2°C cooler than Mg/Ca SST_Gr during BA and are close to U^K′₃₇ SST reconstructions from H1, to Early Holocene, instead.

This indicates that SIMMAX SST reconstructions react in a different way to AF position changes than the other proxies. This could be caused either by an underestimation of Mg/Ca SST_Gr, Mg/Ca SST_Gb and U^K′₃₇ SST amplitudes or by an overestimation of the cooling amplitude by SIMMAX SST.

An overestimation of the alkenone SST has been proposed in previous studies: A similar change in the alkenone signal with respect to SIMMAX SST has been observed in records from the Mediterranean Sea, the eastern African margin, and in the subpolar North Atlantic (Chapman et al., 1996; Sikes and Keigwin, 1996; Sicre et al., 2013), and has been explained with a change in growth period of the alkenone producing coccolithophores from late spring to summer. A change in the coccolithophore growth period could be explained by changes in the water mass properties (e.g., nutrient concentration, temperature, stratification) associated with (1) overall climate forcing such as a change in insolation (Leduc et al., 2010), or (2) a change in the position of hydrographic fronts (i.e., the AF).

(1) Previous work suggests a change in seasonality of the coccolithophores from summer to spring conditions during periods of insolation maxima (Leduc et al., 2010). The time period between H1 and Early Holocene coincides with the period of maximum summer insolation (Berger and Loutre, 1991; Figure 4). To explain the record at the AF, a change of the coccolithophore bloom from spring to summer during the insolation maximum and back to spring conditions after the summer insolation maximum would contradict this concept and seems unlikely. Furthermore, the Holocene cooling trend that is observed in subpolar North Atlantic U^K′₃₇ SST records and has been associated with insolation changes (Leduc et al., 2010) is missing in the Azores record.

(2) Southward displacements of the AF position are indicated by both decreased G. ruber abundances (Figure 3A), and elevated primary productivity (PP) during H1 and YD (Schwab et al., 2012). The pattern of change matches the succession of the deglacial warm and cold events and points to a southward change of the AF during deglacial cold events. The divergence of the SIMMAX SST_su and U^K′₃₇ SST record appears between the Early and Mid-Holocene. This timing does not match with the strong deglacial oscillations in the AF position.

In addition to the major displacements of the AF during cold events, an increased influence of nutrient rich ENACW on the coring site could have had an impact on the phytoplankton bloom and facilitate a late plankton bloom and thus change in the U^K′₃₇ SST seasonality from summer to spring. An increased abundance of N. incompta (10–30%) (Figure 3E) is observed at the coring site from H1 to 10.5 ka BP. N. incompta typically occurs in high numbers in the modern ENACW (Schiebel et al., 2002a,b). Therefore, it can be assumed that more nutrient rich transitional ENACW prevailed at the coring site during this period coincident with the time period of similar SIMMAX SST_su and U^K′₃₇ SSTs (Figure 4). This may be an indicator of a change in the timing of phytoplankton production. Under modern conditions, the delay in coccolithophore production of positive vs. negative NAO phases in the subtropical North Atlantic is about 3 weeks (Henson et al., 2009) and relates to 0.5°C warmer temperatures at 50 m water depth. When added to the observed U^K′₃₇ SST amplitude of 0.4°C, the overall amplitude of 0.9°C U^K′₃₇ SST could be assumed and would be in the same range as the amplitude observed in Mg/Ca_Gr. Thus, a dampening of the overall alkenone SST amplitudes by a change in phenology is not unlikely for the time period between H1 and 10.5 ka BP.

Yet, the calculated amplitudes would still be much smaller than the SST difference of 2.2°C as reconstructed from SIMMAX SST_su between Holocene and BA. Therefore, the observed match between U^K′₃₇ SSTs and SIMMAX SST_su between H1 and Early Holocene, and its shift toward a match between Mg/Ca SST_Gr and SIMMAX SST_su in the Early Holocene cannot be solely explained by changes of the growth period of the coccolithophores, G. ruber, and G. bulloides. If this would be true, the SIMMAX SST_su record would produce too low SSTs. This could be caused by missing analog situations within the calibration dataset (Storz et al., 2009; Telford et al., 2013) possibly because mixed planktic foraminifer assemblages that represent hydrographic fronts are underrepresented in the SIMMAX calibration dataset (Storz et al., 2009). This becomes evident when regarding the similarity index below 0.9% calculated by SIMMAX for the interval between H1 and Early Holocene. To test the influence of the calibration dataset, SIMMAX calculations additionally were carried out using the updated and extended calibration dataset for the North Atlantic of Siccha and Kucera (2017). The new results are similar to the previous results, indicating that the underrepresentation of hydrographic front environments is not the controlling factor for the dissimilarities. Non-analog situations may occur over the entire down core record from the Azores Front-Current System, due to the inclusion of subpolar elements, and (cold) subsurface elements. Relatively cold SIMMAX SST_su are recorded during periods with elevated relative abundances of N. incompta (10–30%) (Figure 5) and decreased relative abundance of G. ruber. Neogloboquadrina incompta typically occurs in high numbers in the ENACW, and it can be assumed that ENACW prevailed at the coring site during BA. This water mass signal seems to overprint the SST signal in the SIMMAX calculations and leads to an underestimation of the SST at the coring site.

The strong effect of water mass properties other than temperature on the SIMMAX record is tested in the following by comparing the different SST reconstructions during cold periods with even stronger changes in the planktic foraminifer assemblage composition.

5.3. Differences in reconstructed temperatures during cold events

Centennial to decadal-scale cold events such as H1 and IACP are associated with extreme cooling in the SIMMAX SST record (Figure 2A). Especially during the short-term IACP cold event, only the SIMMAX SST_su record shows significant cooling with temperatures 5–6°C colder than U^K′₃₇ SST reconstructions (Figures 2, 4). Unfortunately, low abundances in G. ruber and an overprint by Fe-rich crusts on G. bulloides shells impede the comparison with Mg/Ca SST data during IACP and H1. The mismatch between SIMMAX SST_su and U^K′₃₇ SST might either be caused by a bias of U^K′₃₇ toward warm temperatures, admixture of transported alkenones from warmer regions (Mollenhauer et al., 2005), or non-analog conditions that result in anomalously low SIMMAX SST for the cold IACP event.

An underestimation of the cooling in the U^K′₃₇ SST in comparison to SIMMAX SST_su caused by an extreme increase in the abundance of the 37:4 component during cold events that is not considered in the U^K′₃₇ calculations has been excluded by Schwab et al. (2012) as already discussed above (see section “5.1. Reliability of records”). The weaker cooling suggested by the U^K′₃₇ SST could also be caused by lower sampling resolution compared to that of the foraminifer counts not resolving decadal-scale variability. We tested this hypothesis using the same sampling resolution for SIMMAX SST and the U^K′₃₇ SST record (Supplementary Figure 2). At lower resolution, the IACP cooling in the SIMMAX SST remains evident in two data points, indicating that the missing cooling signal is not caused by resolution.

It has been widely discussed that the alkenone temperature signal is smoothing-out short-term SST variability (Bard, 2001b; Ohkouchi et al., 2002; Mollenhauer et al., 2005; Leduc et al., 2010; Laepple and Huybers, 2014). This smoothing might either be related to the large amount of coccolithophore specimens contributing to the alkenones preserved in the sediments (Leduc et al., 2010; Laepple and Huybers, 2014), while about 400 foraminifers per sample are enumerated for the SIMMAX estimation. The smoothing may also be related to differential transport of coccolithophores and foraminifers with ocean currents (Ohkouchi et al., 2002; Mollenhauer et al., 2005), differential sedimentation, and bioturbation (Bard, 2001b). Parallel changes in coccolithophore (Schwab et al., 2012) and planktic foraminifer assemblages show pronounced abundance peaks of transitional to subpolar species during IACP, and indicate that transport-related biases between coccolithophores and planktic foraminifer are unlikely. The coinciding records as well exclude bioturbation as the cause of differential U^K′₃₇ SST and SIMMAX records. A third explanation arises from the modern structure of the AF. Under modern conditions, only a weak temperature gradient is observed in the surface waters at the AF (Alves et al., 2002). Thus, a strong temperature imprint just might be missing in the Alkenone record. In contrast, the SIMMAX_su SST react very sensitive to changes in the foraminifer assemblage composition that are associated with the southward displacement of the AF.

In the following, we interpret the missing cooling signal of the short-term events in the alkenone record as smoothing and lack of an SST gradient across the front. Still, the probability of an underestimation of SIMMAX SSTs during the cold H1 and IACP event needs to be explained.

The possibility of receiving too cold temperatures with the SIMMAX method due to an underrepresentation of mixed planktic foraminifer assemblages from hydrographic fronts in SIMMAX calibration datasets is likely, as discussed above. These non-analog situations of the down core data from the Azores-Front Current-System are evident from low similarities of 0.7–0.85 of the SIMMAX SST calculations (Supplementary Table 2) over the IACP interval, and need to be regarded with care (Pflaumann et al., 1996). More detailed analyses of the faunal assemblage during this IACP interval may help to understand the cause for the dissimilarities.

The cooling of the SIMMAX SST record (Figure 2A) of the IACP event is paralleled by high abundances of T. quinqueloba (Figure 3F) and the temperate subsurface species G. scitula (Figure 2D; Schiebel et al., 2002a). It is likely that the abundances of these species cause cold SIMMAX SST excursions during the IACP. High relative abundances of T. quinqueloba in the modern ocean with up to 20% are found in the NA (Schiebel and Hemleben, 2000). Relative T. quinqueloba abundances of 30–35% are observed in bottom sediments and sediment trap samples in the Nordic Seas (65°N) and southeast of Greenland (52°N) (Pflaumann et al., 1996; Jonkers et al., 2013). Maximum relative G. scitula frequencies of 7% prevail in sediment samples between 40 and 43°N (Pflaumann et al., 1996). In the subtropical North Atlantic, typical relative abundances of T. quinqueloba are in the range of 10–15% (Schiebel, 2002). High relative abundances of T. quinqueloba either reflect colder conditions or a change in the sedimentation of tests of T. quinqueloba during IACP event. A change in sedimentation of tests of T. quinqueloba could have been caused by changes in the water mass stratification that allowed faster sedimentation of the relatively light T. quinqueloba tests (Schiebel and Hemleben, 2017). These changes may have been related to changes in the AF. Interestingly, SIMMAX reconstructions from the Mediterranean Sea (Sicre et al., 2013) show cold SST_su estimates similar to our coring site, also being associated to abundance peaks of T. quinqueloba.

Enhanced abundances G. scitula indicate a southward displacement of the AF from interglacial to glacial conditions and suggest that the abundance changes were caused by changes in the latitudinal position of the frontal systems of the subtropical NA (Schiebel et al., 2002a). Consequently, oscillations of the AF between more southward and northward positions could explain the frequent changes in SIMMAX SST.

Despite the likelihood that the extreme cooling in the SIMMAX record was caused by a strong imprint of the assemblage change, in the following the plausibility of the amplitude of SST change is tested. For this purpose, we use the assumption that the northern limit of the STG was displaced southward during the deglacial cold H1 and IACP events (Repschläger et al., 2015). Further, we assume that the high relative abundances of G. scitula during IACP were associated with the presence of temperate to subpolar waters at the coring site. The latter might even be coupled with an extreme southward displacement of the Arctic Front (e.g., Bard et al., 1987; Eynaud et al., 2009, and citations therein). These events were most probably caused by massive meltwater outbursts from the Laurentian ice sheets that penetrated the NA (e.g., Obbink et al., 2010).

Assuming that subpolar conditions prevailed at the Azores coring site during cold events, coldest possible SSTs would match SSTs observed in the northern NA. Mg/Ca SSTs derived from G. bulloides and N. pachyderma from the subpolar NA (52 to 62°N) indicate a cooling by 7–8°C from 11–15°C down to 6–7°C for the IACP and YD (Came et al., 2007; Peck et al., 2008; Benway et al., 2010; Thornalley et al., 2010) and SIMMAX SST reconstructions indicate that summer SSTs cooled down by 7°C from 15°C to 8°C at 50°N (site BOFS 5K) after BA (Maslin et al., 1995). During H1, it was even colder with SIMMAX-derived SST_su of 3°C at 50°N (Maslin et al., 1995), 7°C according to MAT, 9–10°C according to U^K′₃₇ SST (Pérez-Folgado et al., 2003) and 10°C (MAT) at 40°N (Chapman et al., 2000). Within all of these records, SIMMAX SST_su indicate the strongest cooling and might be overprinted by water mass changes due to frontal movements as supposed at the AF.

All of these reconstructions show that during H1 and IACP SSTs in the subpolar North Atlantic are even lower than our SIMMAX SST_su of 10°C close to the AF. A southward extension of this water mass could be a potential source for the extreme cooling signals.

In order to comprehensively address the potential cooling signal, we also test the feasibility of the SST estimates for sub-millennial scale cold events by comparing our result with those of climate simulations. In several coupled ocean atmosphere model experiments with a freshwater induced AMOC collapse, SSTs cool down by up to 13°C in the subpolar region, whereas the cooling outside of the subpolar region does not exceed 5°C (Stouffer et al., 2006; Liu et al., 2009; Kageyama et al., 2010). Within the uncertainty of the models a 4–6°C cooling at the Azores coring site is reasonable, but seems to be highly dependent on the position of the AF. This highlights the importance of a reliable and detailed reconstruction in paleoceanography and ocean circulation models.

Similar absolute temperatures of 10–12°C are observed in the SIMMAX reconstructions (Figure 2A) during the two major cold events (H1 and IACP). This scenario is unlikely because of the different climatic background conditions before, during, and after the deglacial warming and ice sheet melting. SST signals of this amplitude are neither observed in the subpolar nor polar NA. This cooling is also not supported by the other SST reconstructions at the Azores coring site. Therefore, we assume that temperature amplitudes of up to 12°C as estimated for the IACP cooling when using SIMMAX SST reconstructions at the Azores coring site are overestimated.

5.4. Younger dryas subsurface influence on the SIMMAX SST record

The SIMMAX, U^K′₃₇ SST, and Mg/Ca SST_Gb signals are decoupled during the YD (Figure 2D). SIMMAX SSTs suggest a short-term warming event shortly after the beginning of the YD cooling, while U^K′₃₇ SSTs and Mg/Ca SST_Gb decrease down to 15°C at the onset of the YD, and remain at 15–16°C (Figure 1C) similar to Mg/Ca T_Gtrunc subsurface temperatures. Thus, either U^K′₃₇ reconstructions do not correctly account for this relatively short warming interval, or YD-cooling is masked in the SIMMAX SST record. SIMMAX SSTs may also reflect subsurface SSTs rather than sea surface conditions during cold excursions due to high relative abundances of subsurface dwelling foraminifers as suggested by Telford et al. (2013).

A strong impact of the temperate surface-dwelling species N. incompta (Schiebel and Hemleben, 2000; Schiebel et al., 2002a; Rebotim et al., 2017) on the SIMMAX SST_su record becomes most evident in the close match between the SIMMAX SST record and the inverse record of N. incompta abundance (Figure 5). This record shows an increasing abundance of N. incompta during the cold YD, whereas the subtropical G. ruber become rare, <4% (Figure 3A). Thus, N. incompta becomes one of the dominant species during the YD and its abundance pattern becomes similar to the pattern of the SIMMAX SST record (Figure 5). During this time-period, stratification of the surface water column was different from the modern, and a stronger increase of subsurface than shallow-dwelling species probably led to an increased effect of subsurface conditions on the SST reconstruction. Thus, SIMMAX SST cannot exclusively be regarded as a surface water signal in this part of the YD record given the complex hydrological situation at the AF. As soon as the G. ruber population recovers in the MH, SIMMAX SSTs seem to reflect surface temperatures again.

U^K′₃₇ SST similar to Mg/Ca SST_Gb during YD, again can be explained by similar seasonality and dwelling depths of the alkenone producing coccolithophores and G. bulloides. Similar Mg/Ca SST_Gb, U^K′₃₇ SST and Mg/Ca SST_Gtrunc may relate to subsurface warming (Repschläger et al., 2015), to a collapse of the winter stratification or to yet unknown changes in the life cycle of G. truncatulinoides.

5.5. Comparison of the minor cooling during EH “flickering” phase with deglacial cold events

A series of minor cold-warm oscillations are evident within the SIMMAX SST record during the EH (Figure 2A). These results are not reproduced when using the combined Atlantic and Mediterranean dataset and indicate, that the “flickering” is caused by non- analog situations in the calibration dataset. Yet, the cold events within these oscillations cannot be related to an increase in frequency of the cold-water species G. scitula and T. quinqueloba (Figures 3C, F), as during H1 and IACP. This “flickering” temperature signal in the SIMMAX SST record in the EH might result from minor changes in the faunal assemblage that are related to changes in the position and structure of the AF. These changes are a combination of strongly varying G. bulloides and G. glutinata abundances, and slowly increasing G. ruber and decreasing N. incompta abundances. Due to underrepresentation of AF assemblages (Storz et al., 2009) in the calibration dataset, these minor compositional changes might lead to a pronounced switch of the SIMMAX SST based either on subpolar or transitional samples on the SST reconstruction. Interestingly, the cold events match the periods of several subordinate meltwater pulses of the Laurentian Ice sheet draining the St. Lawrence River (Teller et al., 2002; Jennings et al., 2015). Therefore, these waters may have had a less strong impact on the STG than the main meltwater drainage events during H1, IACP, and 8.2-ka-BP-Event (Repschläger et al., 2015) that are related to meltwater pulses from Hudson Bay (Teller et al., 2002). This assumption still needs to be tested with independent high-resolution δ¹⁸O_w reconstructions by combined Mg/Ca SST and δ¹⁸O measurements. The change in the G. bulloides vs. G. glutinata abundance might be associated with a change in the seasonal succession of species, potentially associated with a higher tolerance of G. glutinata to lower SST. Alternatively, G. bulloides might lack one of the main abundance peaks in autumn, but which more typically occurs further north (Jonkers and Kučera, 2015).

5.6. Effect of SST on derived local δ¹⁸O_sw reconstructions

In the following, we assess the effect of the combination of temperature reconstructions with δ¹⁸O measurements from different proxies to assess δ¹⁸O_w–ice changes at our core site. δ¹⁸O_w–ice was calculated from the temperature reconstructions from Mg/Ca SST_Gr, Mg/Ca SST_Gb, U^K′₃₇ SST and SIMMAX SST_su, in combination with δ¹⁸O records from surface dwelling G. bulloides and G. ruber (Figure 6). Given the offsets in SST reconstructions revealed by our multiproxy study, the combination of temperature reconstructions from different proxy-bearing organisms with δ¹⁸O measurements can have a severe effect on calculated δ¹⁸O_w–ice records. In general, SIMMAX based δ¹⁸O_w–ice records are more variable than other records and reflect the variability of SIMMAX SST records. The G. ruber δ¹⁸O and Mg/Ca SST_Gr based δ¹⁸O_w–ice reconstructions are similar to those based on SIMMAX SST_su during the Holocene (Figure 6A), while the G. ruber δ¹⁸O and U^K′₃₇ based records are systematically lower in δ¹⁸O_w–ice values. For the YD and BA records based on G. bulloides δ¹⁸O, Mg/Ca SST_Gr, Mg/Ca SST_Gb, and U^K′₃₇ (Figure 6B) are in good agreement with exception of short-term freshening events coinciding with the extreme cooling events in the SIMMAX record.

The average differences between the U^K′₃₇ and SIMMAX based reconstructions amount to 0.5 ‰ SMOW in the Holocene. This difference can be regarded as an estimate for the imprint that seasonality can cause on the δ¹⁸O_w–ice signal. Results of SIMMAX based reconstructions indicate a strong freshening of 2 ‰ SMOW during cold H1 and 8.2-ka-BP-Event, whereas freshening amounts to 1 ‰ SMOW in all other reconstructions. Hence, the too cold SIMMAX SSTs during these intervals lead to an overestimation of the freshening by 1 ‰ SMOW. These results show that differences in the productive seasons of the different organisms carrying the SST or δ¹⁸O signals may result in differences in calculated δ¹⁸O_w–ice at the same order of magnitude as a signal of seawater freshening itself. Therefore, it is important to combine δ¹⁸O signals and SST signals for δ¹⁸O_w–ice reconstructions representative of the same season and dwelling depth, ideally together with Mg/Ca temperatures and δ¹⁸O reconstructions from the same planktic foraminifer sample.

5.7. Reconstruction of seasonal differences using a multiproxy approach

Changes in the phenology of coccolithophores and planktic foraminifers in the North Atlantic are a matter of changes in the growth period (Hopkins et al., 2015; Kretschmer et al., 2018). However, little is known about seasonal changes, i.e., onset and duration of seasons on longer time scales such as glacial-interglacial cycles, although seasonality might have a major impact on and might be impacted by climate variability. Under the assumption that different species of planktic foraminifers and coccolithophores produce Mg/Ca and alkenones signals during different seasons in the North Atlantic, these differences could be used in a combined approach to reconstruct temperature changes between spring, summer, and autumn. The difference between average Mg/Ca SST_Gr (summer SST) and U^K′₃₇ (spring SST) increases from 2.3°C during BA to 2.7°C and 2.9°C over the Early and Late Holocene, respectively. The latter increase (BA to Late Holocene) in SST difference may relate to an increase in absolute seasonal differences in SST between BA and Late Holocene. In addition, extremely cold winter temperatures are proposed for the YD interval (Denton et al., 2005). Colder than modern winters might also have prevailed during BA and may have led to a shift of the growth periods of different species to different degrees. Consequently, the spring bloom may have been delayed while summer conditions were possibly constant during peak summer SST.

6. Conclusion

Distinct differences in multiproxy SST records are caused by the unique nature of the different archives and are influenced by the ecology of the proxy bearing organisms, including depth habitat and main growing season. Once understood, these differences may help to better reconstruct changes in seasonal temperature and surface ocean stratification, which are major unknowns in global climate reconstruction.

In order to better understand the different imprints on the SST proxies in the North Atlantic, a multiproxy SST record from the boundary between the subtropical and the temperate North Atlantic was analyzed and discussed under consideration of the phenology and habitat depth of different paleoceanographic proxies. The comparison between the SST indicates that the different temperature proxies of planktic foraminifer based transfer calculations (SIMMAX), Mg/Ca of different planktic foraminifer species, and alkenones do not only record changes in annual temperatures but also give an insight into changes in different seasonal temperatures. The main findings are that

• Mg/Ca SST reconstructions from G. ruber are similar in the BA and Holocene, and most likely reflect summer SSTs down to 50 m water depth at the northern edge of the subtropical gyre.

• Mg/Ca SST reconstructions from G. bulloides mainly reflect temperatures of the upper 60 to 100 m water depth in spring, and, thereby, in the same season as the alkenones.

• Mg/Ca T reconstructions from G. truncatulinoides reflect temperatures at thermocline depth (100–150 m) that are decoupled from the sea surface and indicate subsurface warming during BA ad YD period.

• U^K′₃₇ SST reflect surface water temperatures in the upper 50 m during late spring and show the well-known deglacial North Atlantic succession of warm and cold events. U^K′₃₇ SST may underestimate the overall deglacial temperature amplitudes and the magnitude of short-term cooling events (e.g., H1, IACP, and 8.2-ka-BP-Event) that interrupt the long-term deglacial and Holocene SST evolution due to changes in phenology.

• Faunal assemblage changes seem to record water mass changes, even in situations of minor SST changes, which are not unequivocally displayed by SIMMAX SST reconstructions. Thus, SIMMAX SSTs using the original calibration dataset (Pflaumann et al., 1996) are more sensitive to changes in water mass distribution rather than changes in the absolute SST, and result in over- or underestimation of SST changes. These differences become greater at strong SIMMAX SST cooling during the two major cold events H1 and IACP, and minor cold oscillations during the Early Holocene. An overestimation of SSTs during the early YD is probably caused by SIMMAX SSTs reflecting subsurface conditions that indicate thermocline warming.

• Despite methodological differences, strong cooling at the core site SW of the Azores during the major cold events H1 and IACP is indicated by the occurrence of subpolar faunal elements.

• Differences between U^K′₃₇ SST and Mg/Ca SST_Gr reconstructions could be used to assess changes in duration and/or amplitude of seasonal changes in our study region.

• Ecological differences resulting from seasonality and dwelling depth in the chosen SST signals can have strong effects on δ¹⁸O_w–ice reconstructions, probably of an amount similar to the effects of freshening events.

Data availability statement

The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

JR, RalphS, and MW contributed to conception and design of the study. JR and MW performed the faunal assemblage counting and SIMMAX calculations. JR analyzed the data wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

Publication of this manuscript has been supported by the European Science Foundation (ESF) under the EUROCORES Program EuroMARC, ERASCT-2003-980409 of the European Commission, DG Research, FPG. Financial support was provided by the German Research Foundation (DFG grant number WE 2679/6-1).

Acknowledgments

The results represented here are obtained from cores taken on board the R/V Marion Dufresne (French Polar Institute, IPEV) during the IMAGES-AMOCINT MD168-cruise in the framework of the 06-EuroMARC-FP-006 Project. We thank the captain and crew, as well as our colleagues on board for their help. Furthermore, we also thank Bernard Dennielou and IFREMER for providing us with material from core GEOFAR KF16.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1176278/full#supplementary-material

References